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Non-Radiative Test of a New SPH Scheme 
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ABSTRACT 

We adapt new schemes of gravitational N-body dynamics and smoothed particle hy- 
drodynamics (SPH) to the tree N-body/SPH galactic chemodynamics code GCD+. The 
applied new schemes include the adaptive softening suggested by Price & Monaghan 
(2007), to improve the self-gravity calculation, and artificial viscosity and thermal con- 
ductivity suggested by Rosswog & Price (2007) and Price (2008), to model discontinu- 
ities and Kelvin- Helmholtz instabilities more accurately. We first present non-radiative 
test simulations and contrast the results with those obtained with conventional hydro- 
dynamical schemes. The new schemes lead to a significant improvement in the ability 
of TreeSPH codes, such as GCD+, to capture strong shocks and model Kelvin-Helmholtz 
instabilities. In light of the well-known discrepancy between entropy profiles predicted 
by SPH codes and those from adaptive mesh refinement (AMR) codes, we next ex- 
plore the impact of these new hydrodynamical schemes upon cosmological simulations 
of galaxy cluster formation. We demonstrate that our new version of GCD+ leads to 
higher entropy in the cluster centre, with an associated entropy profile consistent with 
that obtained with AMR codes. We suggest that the newly-applied SPH scheme may 
ameliorate this classical discrepancy between SPH and AMR codes. 
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1 INTRODUCTION 



Since it was int roduced 
iGingold &; Mon aghan 



by iLucvl (|l977h and 
the smoothed particle 
hydrodynamics (SPH) methodology has become a reg- 
ular tool for the numerical si mulation of a wid e range 
of astronomical phenomena. iHernquist &; Katzl (p.989) 
were the first to suggest that the SPH approach would 
also prove invaluable in the simulation of galaxy for- 
mation and evolution. Since then, a number of SPH 
codes have been developed to simulate such systems, 
incorporating various physical processes ranging from 
radiative cooling to star formation and supernovae 
(SNe) feedback /e.g. iKatd Il992l: iNavarro fe Whit d [l 993: 



1999 



2004 



Katz et al. 1996; S teirimetz fc M uller 1995; M ori et all 1999 
Carr aro et al.l 1 19981 : iKawatal Il999l : ISommer-Larsen et al 



Springel et al.ll200ll;lKobavashill2004l;lGovernato et al 
Okamoto et alJl2008l : ISaitoh et alJl2008h 



Parallel to the development of such particle-based 
codes, grid- or mesh-based approaches have been employed 
for modeling the formation and evolution of galaxies (e.g. 
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ICen fe Ostrikeilll992h . Algorithmic enhancements to a fixed 
grid approach, such as adaptive mesh refinement (AMR), 
has led to a massive improvement in the capability of grid- 
based codes for simulations which require a large dynamic 



range, includin g those of galaxy formation (iTevssierl 



Kravtsovl [20031: iTasker fe Brvanl [20081 : IGibson et al. 



2002 



2008 



Joung et alj I2008T). Code comparisons between SPH and 
AMR (j[r enk et al. 199jl lAscasibar et al.l 120031 : IVoit et al 



| 2005| : lO'Shea et al l| 2005 [ : IGibson et aT1 l2008: Mit chell et al 
2008: ITasker et al.ll2008t ) demonstrate that the competing 
approaches lead to generally consistent results. That said, 
comparing the results of non-radia tive s imulations of the 
formation of a galaxy cluster, iFrenk et al.l ([19991 ) claim that 
SPH codes lead to lower e ntropy in the central region of the 
simu l ated cluster (see also lAscas ibar et al. 20 031: IVoit et al. 
| 2005| : lO'Shea et al.l|2005l: iDolag et al] 120051 : iMitchell et al. 
booalWadslev et alJfeOOa ). They suggest that SPH may un- 
derestimate turbulence in the central region. 

lAgertz et all (|2007h carried out a series of experiments 
in order to compare and contrast SPH and AMR in more 
of a "controlled" environment. They conclude that there 
is a "fundamental" discrepancy between these approaches, 
by demonstrating that SPH, at least in its conventional 
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form, cannot capture Kelvin-Helmholtz instabilities (KHI) 
as accurately as an AMR approach. They further suggest 
that this discrepancy is not due to one of resolution, but 
is a fundamental attribute of the scheme its elf (see also 
llmaeda fc Inutsukalliool lOkamoto et al.ll2003h . 

Recently, IPricd (|2008l ) described a new scheme to im- 
prove the conventional implementation of SPH and demon- 
strated the successful capture of KHI. In what follow, we 
apply this schem e to our origina l galactic chemodynam- 
ics code, GCD+ (|Kawata & Gibson! [2003a ). GCD+ is a three- 
dimensional tree iV-body/SPH code that incorporates self- 
gravity, hydrodynamics, radiative cooling, star formation, 
supernova feedback, and metal enrichment. At its heart, the 
new scheme differentiates itself from the conventional ap- 
proach via the manner by which diffusion of thermal energy 
is introduced; we adopt pr imarily the formalism described 
bv lRosswog &; Pried (|2007l ). We demonstrate that this new 
scheme does indeed advance the abilities of SPH codes in a 
suite of controlled tests (some of which have been described 
by Price 2008), in addition to leading to more realistic sim- 
ulations of galaxy formation and evolution. To our knowl- 
edge, our work is the first demonstration of this new scheme 
within the context of cosmological simulations which follow 
both gas and collisionless dark matter (DM) dynamics simul- 
taneously. We focus only upon non-radiative simulations in 
this paper, and will study cases including radiative cooling 
and star formation in a forthcoming work. 

Section [2] describes briefly the implementation of this 
new scheme within GCD+, hereafter differentiating this ver- 
sion from the conventional one by the moniker "GCD+ Mkll" . 
Section [3] compares the performance of GCD+ Mkll with the 
conventional version of GCD+ under a range of test conditions, 
summarising our results in Section [4] 



2 GCD+ MKII: ADVANCING GALACTIC 
CHEMODYNAMICS 

We now describe the specific modifications made to the 
galactic chemodynamics code GCD+, which themselves are 
patterned closely after the methodology described by 
IRosswog fe Price! (|2007h . As such, we only outline the fi- 



nal formulae adopted, and refer the interested reader to 
IRosswog fe Pried |2007) for their formal derivation. 
The density of the i-th SPH particle is defined by 



y^mjWjnj, hi 



(i) 



w here W i s the S PH smoothing kernel seen in equation JT} 
of iKawatal ([l999), nj = \xi — Xj\, and hi is the smooth- 
ing length of the i-th particle. We note in passing that 
GCD+ Mkll only takes into account the smoothing length of 
the i-th particle, hi, to derive the density, while the origi- 
nal version of GCD+ ([Kawatal 1 19991 ) used the pair-averaged 
smoothing length, hij = (hi + hj)/2. The smoothing length 
is determined by 



hi 



1/3 



(2) 



where we adopt 77 = 2.4. Note that in our definition of the 
kernel, our smoothing length corresponds to twice that used 



bv IRosswog fc Pried <|2007h , who adopt rj = 1.2. The solu- 
tion of equation J2} is calculated iteratively until the rela- 
tive change between two iterations is smaller than 10 _ 3 (see 
IPrice &; Monaghanll2007l for more details). Euler's equation 
is written as 



dt 



/ ^v,u-, J (/ ) ,)--^v,u ./, 

Pi lL jpj 



+ Qv„ 



- G^rriji-^ — J - >e tJ 

- ^J2 rr ^ {^ iWij{hi) + ^ ViWij{hj) ] ■ (3) 

The first term of equation J3} corresponds to the pres- 
sure gradient, where Wij(hi) = W(rij, hi), ViWij(hi) = 
dW(rij,hi)/dxi and 



dhi dWik(hi) 



dpi ' dhi 



(4) 



From equation (|2]), dhi /dpi = —hi/ (3 pi). The second term 
of equation J3]) corresponds to the artificial viscosity ( AV) , 



Qv,i 



Vij • eij 



Pij 



ViWij (if Xij • Vi 



= (otherwise), 



where Vij = V{ - Vj, eij — (xi — Xj)/\xi - Xj\, pij — (pi + 
pj)/2 and 



2 \Qi 



The signal velocity v s i g adopted is 

Cs,i ~h Cs,j SVij • Sij 



(6) 



(7) 



where c s ,i is the sound velocity of the i-th particle. The 
amount of AV is controlled by a time-dependent parameter, 



at/ = ^(at v (t) + afmfi + f 3 ), (8) 
where 

f = Kv^ (9) 

l(V-v)»| + |(V x v)i\ +0.0002c Sji /fe' v ; 

in order to suppress AV in pure shear flows. The viscous 
parameter a^ v (t) varies with time. IRosswog fc Pried (|2007h 
suggested the foll owing functi on to evolve this viscous pa- 
rameter (see also lMorrislll997l ): 



AV AV 

ot i a i — ao 



dt n 

where ao — 0.1 and 

_ hi 
n ~ 0.2c sA ' 



+ Si, 



IRosswog fc Pried (|2007h adopt the source term, 
Si = max(-Vi • «i,0)(2 - af v ), 



(10) 



(11) 



(12) 



and set the maximum of the viscous parameter to be a max — 



Non-Radiative Test of a New SPH Scheme 3 



2. We call the version of our code which adopts this partic- 
ular source term "Mkll+lowa", hereafter. As shown later, 
we found that this formulation led to insufficient viscosity 
for the capture of some strong shocks. As such, our fiducial 
model uses 



Si 



(V,-^) 2 (4-«f v ) (if Vr^ <0) 



(otherwise) 



(13) 



and employs a„ 



4. We call this fiducial code 



"Mkll+higha" . The third term of equation © corresponds 
to the gravitational force, and employs the adaptive grav- 
itation al force softening suggested in IPrice &; Monaghanl 
(|2007T ). where the softening length is matched to that of 
the smoothing length. The fourth term of equation (|3} is 
the correction term for adaptive softening, where 



r _ dhi sr^ d<j)ij{hi) 



(14) 



We apply a cub ic splice softening, as suggested by 
|PrjcejLMonagh ar ] (|2QQ7h : the associated formulae for <b 
and d<j)/dh can also be found in their paper. 
The energy equation is written as 



dui 
~dt 



3 



rrijVij • ViWijihi) 



+ Qk 



(15) 



where m is the thermal energy of the i-th particle, and Q u ,ij 
is zero if Xij • Vij > 0. Otherwise, it is described by 



Qu,i — 



C i 
Ctij (Ui 



(16) 



where afj = (af+af)/2. The second term within the paren- 
theses of equation (]T6l) corresponds to the artificial th ermal 
conductivity (AC) jllosswog fe Pricell2007l : |Pricdl2008h . The 
thermal conductivity parameter, a c , evolves following 

dou 



at n 
where the source term is 



(17) 



of higher mass. Note that for the DM we use the mass den- 
sity, pi,r>M = Sj??ij,DMM / (rij, hi) j in equation J2]). Therefore, 
this higher ?7dm results in a DM softening length which is 
larger than that for the gas when the number density of DM 
particles is comparable to that of the gas. We also set the 
minimum softening length e m i n for both N-body and SPH 
particles to be e m i n = (c e ra p (MQ)) 1//3 pc, where m p is the 
mass of the particles. For the work described here, we em- 
ploy c e — 3000, which results in similar softening to t hat 
we used in our previous work (jKawata fe Gib son 2003a .b). 
When hi — e m ;n, we set Vt t — 1 and = 0. 

For SPH particles, we apply an individual timestep 
scheme to integrate equations J3} and JT5J, and for col- 
lisionless particles we use equation (l2 0]l. We also em- 
ploy the timestep limiter suggested by ISaitoh &; Makinol 
(2008). The timestep for SPH particles is based upon 
dU = mm(dtcFL,ij,dtr>YN,i)i where the Courant-Friedrich- 
Levy condition is calculated by 



dt c 



FL,ij 



0.5/i 2 

Vdt,ij 



(21) 



where vat,ij 



^sig, ij if Xij 



Vij < 0, otherwise Vdt,ij = 
^ • eij). We set Ccfl = 0.2 for Mkll+alow. 
On the other hand, Ccfl = 0.1 is adopted for Mkll+ahigh, 
as we found a lower Ccfl was required in a one-dimensional 
shock-tube test. The requirement that the force should not 
change significantly within one timestep is satisfied by 



dtr>YN,i = Cdyn 



O.bhj 
\dvi/dt\ 



1/2 



(22) 



We set Cdyn = 0.2 for SPH particles. The timestep for 
collisionless particles is derived solely from equation ([22]) . 
and we adopt Cdyn = 0.2. 



3 RESULTS 

Having outlined the improvements made to GCD+, we now 
compare its performance with that of the original version. 
As noted in Section [2] we explore the impact of two differ- 
ent forms for the AV source term: Mkll+lowa (using equa- 
tion [12]) and Mkll+higha (using equation 1 13[) . 



and 



omhiWui 



V Ui — 2 > vfij — 



pj 



(18) 



(19) 



For collisionless N-body particles (such as DM), only 
the gravitational force term of equation © is taken into 
account, as follows: 



dvj 
dt 



-G^m 3 



3 



^ij{hi) +(/>ij(hj) 



^ViWijihi) + ^ViWijihj) } .(20) 



The adaptive softening length of the DM is calculated us- 
ing equation J2]) and for the DM we apply ?7dm = 4, as the 
simulations under consideration generally have DM particles 



3.1 One-dimensional shock-tube test 

Not surprisingly, our first experiments involve a version of 
the classical one-dimensional shock-tube test. The initial 
conditions are set by assuming the simulation region spans 
from x — —0.5 to x = 0.5; the region for which x < is set to 
(pi, Pi, vi) = (1,1,0), using 960 particles, and the region for 
which x > is set to (p 2 ,P2,v 2 ) = (0.125,0.1,0), using 120 
particles, adopting a 7 = 5/3 throughout. Figs. [T] [2j and [3] 
show the results for the original version of GCD+, Mkll+lowa, 
and Mkll+higha, respectively. One can see a clear jump in 
thermal energy and pressure at contact discontinuities in the 
original version, while in the Mkll+lowa version, the con- 
tact discontinuity is resolved, and the pressure and thermal 
energy distribution is much smoother. It is also remarkable 
that the number of particles employed to resolve the shock 
front in this case is so low. Although less smooth compared 
with that of the Mkll+lowa version, the Mkll+higha im- 
plementation also resolves the contact discontinuity better 
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0.5 -0.5 



Figure 1. Results of a one-dimensional shock-tube test with the 
original version of GCD+ at t = 0.2. The grey line represents the 
analytic solution. 




0.5 -0.5 



Figure 2. As in Fig. ^ but now using the Mkll+lowai version 
of GCD+. 
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Figure 4. Radial density (upper) and pressure (lower) distribu- 
tions at t = 0.04 in the point-like explosion test using the original 
version of GCD+. The grey line shows the analytic solution. 



than that found using the original GCD+. However, because 
we assume a higher a AV in this case, a greater number of 
particles is required to resolve the shock front. A smooth 
pressure distribution at the contact discontinuity is the key 
to simulate accurately KHI (Price 2008); as such, it would 
appear that both the Mkll+lowa and Mkll+higha versions 
of GCD+ are promising tools for modeling KHI within an SPH 
framework. 



T 




-0.5 0.5 -0.5 0.5 

x x 



Figure 3. As in Fig. [T] but now using the Mkll+highai version 
of GCD+. 



3.2 Point-like explosion test 

We next consid er the Sedov- Taylor- type sph erical explosion 
test. Following ISoringel fe Hernauistl (2002), we set a peri- 
odic boundary box with low-temperature and homogeneous 
density (p = 1), using 64 3 SPH particles. At t = 0, we de- 
posit E — 1 energy on the central particle and simulate the 
evolution thereafter. The analytic solution can then be de- 
rived via the adoption of Sedov- Taylor self-similarity. Fig. [4] 
shows the result at t — 0.04 using the original version of 
GCD+, while the grey line represents the analytic solution; it 
is clear that the original version of the code was incapable of 
capturing a strong shock such as this. The Mkll+lowa ver- 
sion show much better performance in this regard (Fig. [5}, 
although even still, particles do "penetrate" the shock front. 
Conversely, the Mkll+higha version does capture success- 
fully the associated strong shock, as seen in Fig. [6] As such 
strong shocks can be an important byproduct of the pro- 
cesses governing galaxy formation and evolution, we hence- 
forth adopt the Mkll+higha version as our fiducial code. 
Fig. [7] demonstrates that the energy conservation form of 



Non-Radiative Test of a New SPH Scheme 5 



— i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — r~ 




i i \ i I TT*i"~ i I i i i i I i i i i I i i i i 



0.1 





0.04 



Figure 7. Total energy as a function of time for the original 
version of GCD+ (dotted line) and the new Mkll+higho; version 
(solid line), for the point-like explosion test. 



the new SPH scheme described in Section [2] can conserve 
the total energy even in this strong shock case (see also 
IPrice &; Monaghanll2007h . While not shown, the Mkll+lowa 
version also shows the same level of energy conservation. 



Figure 5. As in Fig. [4] but now using the Mkll+lowo; version 
of GCD+. 
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Figure 6. As in Fig. [4] but now using the Mkll+highai version 
of GCD+. 



3.3 Self-similar collapse test 

IBertschingerl (|l985h derived a self-similar solution for the 
collapse of an overdens e perturbation in an E instein-de Sit- 
ter (Q = 1) Universe. [Navarro fc White! ([1993) introduced a 
tes t simulation based upon this self-similar solution. Follow- 
ing [N^rro^Whiti {l993) , we consider a spherical volume 
which initially follows the Hubble expansion, and set a cen- 
tral spherical perturbation with mass of 0.05M to t and radius 
of 0-l.Rini, where M to t is the total mass and R in i is the initial 
radius of the simu l ation sphere. We use more particles than 
iNavarro fe White! (jl993T ). adopting n h = 15515 for SPH par- 
ticles and riDM = 15515 for DM. We also assume Qb = 0.2, 
and therefore the DM particles are four times more massive 
than the SPH particles. 

The dimensionless parameters for radius, A, radial ve- 
locity, V r , density, D, and pressure, P, are defined by 

K(A) = —vr(r,t), 
D(A) = 

PH 

P (A ) = (1\ 2 P±A. (23 ) 

Figs. [8] [9] and [10] show the results in these dimensionless 
parameters. Dots within each panel represent the simula- 
tion results for the gas particles, while the grey line s corre- 
spond to the analytic solution of iBertschingerl (|l985h . Fig. [8] 
demonstrates that although the original version of GCD+ 
roughly reproduces the analytic solution, the radial velocity 
has too much scatter around the shock front, and a small 
fraction of particles obtain significantly high pressure at a 
radius just outside the shock front (log A ~ —0.35). 

The Mkll+lowa version shows much less scatter in the 
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Log A 

Figure 8. Normalised velocity (upper), density (middle), and 
pressure (lower) distribution at an arbitrary time, employing the 
Bertschinger (1985) test with the original version of GCD+. The 
grey line represents the analytic solution. 

velocity distribution (Fig. [9]), with the shock front signifi- 
cantly more apparent. The Mkll+higha version possesses 
even less scatter in the velocity distribution, due to the 
higher AV, although the "sharpness" of the shock front ap- 
pears comparable to that seen with Mkll+lowa. It would 
appear that in three-dimensional simulations, the high a AV 
scheme may not change the resolution of the shock signifi- 
cantly. 

3.4 KHI test 

lAgertz et all (|2007h introduced a straightforward test which 
allows for a given code (particle- or grid-based) to be as- 
sessed in terms of its ability to resolve KHI (see also IPricd 
2008). In this section, we demonstrate that our fiducial 
model (Mkll+higha) is now succes s ful in resolving such 
instabilities. Following lAgertz et al. I <l2007h . we consider a 
periodic boundary box with x — —0.5, 0.5, y — —0.5, 0.5, 
and z = —1/64,1/64. The region within \y\ < 0.25 is set 




-2 -1.5 -1 -0.5 0.5 

Log A 



Figure 9. As in Fig. [8] but now using the Mkll+lowa version 
of GCD+. 

to be the high-density region (with ph = 9.6), while the 
rest is the low-density region (with p\ — 1). Equal mass 
particles are adopted in both regions, and 544 (256) parti- 
cles are used to cover the x-axis for the high-density (low- 
density) region. The two regions are in pressure equilibrium 
and we assume Ph = P\ = 2.5. The high- density (low- 
density) region has velocity v x = —0.5 (v x — 0.5). We also 
added sinusoidal perturbations on the vertical velocity, us- 
ing v y (x) = 5v y sin(A27rx), setting Sv y = 0.025 and A = 1/6. 
These perturbations are only imposed to the particles within 
\y\ = 0.225,0.275. As before, we assume 7 = 5/3; this con- 
dition leads to a timescale for KHI of tkhi = 0.57 

Fig. [11] shows the density and pressure distribution at 
t — tkhi for the Mkll+ce high version; t his c an be com- 
pared directly with Fig. 13 of lAgertz et al 1 (|2007h . The latter 
demonstrates that instabilities do not appear in their SPH 
simulations. On the other hand, simulations with the grid 
codes employed capture the deve lopment of KHI wit hin the 
expected KHI timescale, leading lAgertz et al.l (|2QQ7h to ar- 
gue th at this is a fundament al problem for conventional SPH 
codes. lOkamoto et al.l (|2QQ3h also show in their Fig. 9 that 
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Figure 12. Density (left) and pressure (right) distribution at 
t = 2tkhi for the blob test, using the Mkll+higho; version of 
GCD+. The initial blob consists of 10 6 particles. 



Figure 10. As in Fig. [8] but now using the Mkll+highai version 
of GCD+. 
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-0.4 -0.2 0.2 



Figure 13. As in Fig. 1121 but for a lower-resolution simulation. 
The initial blob consists of 10 5 particles. 



although SPH can develop weak KHI, they quickly disap- 
pear owing to numerical momentum transfer. Fig. [11] demon- 
strates that the new version of our SPH code is capable of 
capturing KHI, and leads to similar re s ults t o those found 
using the AMR codes in lAgertz et all (120071). We confirm 
that the new SPH scheme proposed by iRosswog &; Price! 
(2007) appears t o solve this apparent fundamental problem 
of SPH (see also IPricd 2008. who demonstrates it with two 
dimensional test). 



Figure 11. Density (left) and pressure (right) distributions at 
t — ^kh within the KHI test, using Mkll+highai version of GCD+. 



3.5 Blob test 

In the previous section, we demonstrated that the new ver- 
sion of our SPH code - i.e. Mkll+higha - is capa ble of cap- 
turin g KHI in a simple shear flow simulation. lAgertz et al.l 
(2007) also introduced a so-called "blob test" , a more appro- 
priate analogue f or simulating galaxy formation and evo- 
lution. Following lAgertz et al.l (|2QQ7h (Oscar Agertz, priv 
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Figure 14. As in Fig. [12] but employing the original version of 
GCD+. 
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Figure 15. As in Fig. 1131 but employing the original version of 
GCD+. 

coram), we set a high-density spherical "blob", with tem- 
perature, Tci = 10 6 K, density, p c i = 0.07281 M kpc -3 , 
and radius 193 kpc, within a low-density field of tempera- 
ture Text = 10 7 K and density /w = 0.007281 M kpc -3 
(note: the two regions are in pressure equilibrium). We con- 
sider a periodic box with {L x , L y , Lz} = {1, 1,3 Mpc}. The 
blob is set at {0, 0, —0.75 Mpc} and initially assumed to be 
static, but the low density field is assumed to have a ve- 
locity of v z = 1000 km s _1 . The particles with the same 
mass are distributed randomly so that the above condition 
is satisfied. 

Fig. [12] shows the results at t ~ 2tkhi for the simula- 
tions with the new Mkll+higha version, adopting 10 6 parti- 
cles within the blob. One can see that the blob is disrupted 
significantly due to the instability. We should note that dis- 
ruption is reduced relative to th e results seen with AMR 
codes (see Fig. 4 of lAgertz et al. I l2007h : as such, it may be 
the case that our new code still underestimates the effects 
of KHI. 



t 1 1 1 1 1 1 1 1 1 1 1 1 r 
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Figure 16. Gas density (upper), temperature (middle), and en- 
tropy (lower) profiles at z = of the SBC simulation, using the 
original version of GCD+ (open triangles), the Mkll+higha; version 
(solid circles), and the Mkll+higha' version without AC (crosses). 



Fig. ITTTI shows the same simulation, but with lower res- 
olution (where the blob consists of 10 5 particles). Compar- 
ing with Fig. 1121 the low-resolution simulation significantly 
underestimat es the effects of K HI. Note that the blob test 
in Fig. 4 of lAgertz et al.1 (|2Q07h has the same resolution as 
our low resolution simulation. We also carried out a simu- 
lation with the sa me initial condition s (kindly provided by 
Oscar Agertz) as lAgertz et alj (|2007h , however the results 
were comparable to their published SPH results. We con- 
clude that even with our new version of GCD+, in order to 
accurately capture KHI, high- resolution is required in both 
the low- and high- density regions. One must be ever- vigilant 
when it comes to the appropriate resolution when simulating 
systems where KHI become important. 

For reference, we also show the results employing high- 
fFig.H4]) and low- (Fig. [15} resolution simulations, but with 
the original version of GCD+. Surprisingly, the original ver- 
sion show significantly more disruption than the new ver- 
sion. However, as shown in Sections 13.21 and 13.31 the original 
version employs an artificial viscosity which is too small, 
which is the agent responsible for allowing the instability to 
become more developed. 
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Figure 17. Gas density (left), temperature (middle), and line-of-sight velocity (right) distributions in an arbitrary slice through the 
centre of the simulated cluster, using the original version of GCD+. 
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Figure 18. As in Fig. [TT1 but employing the Mkll+highai version of GCD+. 
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Figure 19. Gas temperature vs. density map for the original version of GCD+ (left), the Mkll+higho; version (middle), and the 
Mkll+higho; version without AC (right). 
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3.6 Santa Barbara Cluster 

Our final test is based upon the so-called Santa Barbara 
Cluster (SBC) simulation, a now-standard test for cosmo- 
logical hydrodynamics codes (Frenk et al. 1999). The ini- 
tial conditions correspond to a cosmological simulation in 
which a galaxy cluster forms at the centre of the final vol- 
ume. We extract a spherical region with a radius of 32 h~ x 
Mpc which is expanding with the Hubble flow, and assume 
an open boundary. The gas and DM particle masses in the 
highest-resolution region are 8.67 x 10 8 and 7.80 x 10 9 M , 
respectively. We run the simulation using both the origi- 
nal version of GCD+ and the new Mkll+higha version, and 
compare the results. The minimum softening (and there- 
fore smoothing) lengths for gas and DM particles are set to 
13.8 and 28.6 kpc, respectively, for the Mkll+higha version. 
Conversely, the original version of GCD+ uses fixed Plummer 
softening lengths of 11.6 and 23.0 kpc, respectively, for gas 

and DM particle s. 

Frenk et al.l (1999) demonstrate that the various codes 
employed in their analysis lead to generally consistent re- 
sults. However, they also found a clear difference in the 
entropy profile between the particle- and grid-based codes, 
with the latter resulting in systematically higher entropy 
than the SPH codes, where entropy is defined as s = 
ln(T/p 2 / 3 ). We compare density, temperature, and entropy 
profiles at z — in Fi g. ITGl Comparing with Figs. 12, 17, 
and 18 of I Frenk et al.l ([ 19991 ), the original version of GCD+ 
shows c onsistent results w ith those from the SPH codes 
used by iFrenk et al.l |l999). It is very interesting that the 
Mkll+higha version of GCD+ shows significantly different re- 
sults from the original version of GCD+. The Mkll+higha 
version results in a much lower density, higher temperature, 
and hence significantly higher entropy at the centre of the 
cluster. As a result, the profiles of the Mkll+higha version 
of the simulatio n are s imilar to those for the grid codes used 
by IFrenk et "ah I (E9990 

In particular, our results show a 
marked similarity to those of the AMR code labeled therein 
as "Bryan" (a precursor to what is now referred to as Enzo) . 

We also ran a simulation with the Mkll+higha version 
of GCD+, but without AC, the results of which are shown as 
crosses in Fig. [16] Such a test leads to a result similar to that 
seen when using the original version of GCD+. We conclude 
that it is AC which ultimately leads to the difference seen in 
the profiles generated with the original version of GCD+ and 
the new Mkll+higha version. 

Figs. [TT1 and [18] show the density, temperature and 
line-of-sight velocity distributions, in the central 4 Mpc 
region of the cluster. We can see that the original ver- 
sion of GCD+ shows more small-scale perturbations, partic- 
ularly in the temperature distribution. On the other hand, 
the Mkll+higha version results in a much smoother dis- 
tribution. Finally, we compare the temperature and den- 
sity distributions for the gas particles within virial radius 
in Fig. 1191 As expected, the Mkll+higha version shows a 
much smoother distribution. Differences in the temperature 



1 Note that it is still un clear what is th e correct entropy profile 
for the SBC. For example, Springel (2009) demonstrates that their 
moving mesh code is capable of modeling KHI, while the code 
leads to similarly lower entropy at the centre of the SBC to what 
conventional SPH codes predict. 



and density distributions are seen in the simulations of Sec- 
tion [3]3] Therefore, we expect the same results even in the 
case of spherically-symmetric monolithic collapse - i.e., in 
the absence of hierarchical clustering. We conclude by not- 
ing that we have still lack a clear understanding as to why 
the new version of GCD+ results in entropy profiles which 
differ from conventional SPH codes. 



4 SUMMARY 

We implement new AV and AC schemes within an N- 
body/SPH galaxy formation code, applying adaptive soft- 
ening to both the N-body and SPH particles. In this paper, 
we study how these new schemes work within the context 
of non-radiative simulations, and focus on the effect of the 
combination of the new AV and AC. 

We demonstrate that this combination helps to 
"smooth" the thermal energy at the contact discontinu- 
ity. Because of this improvement, the new code succeeds 
in capturing KHI in a manner that conventional SPH treat- 
ments do not. In essen ce, this confirms that th e A V and 
AC sc heme proposed by iRosswog fe Pried |2007) and IPricd 
(2008) remedies the fundamental problem of SPH outlined 
bv lAgertz et al.l (|2007T ) . It is also shown that to properly sim- 
ulate the KHI, resolution remains critical, even for our new 
implementation. Therefore, one must remain vigilant when 
it comes to resolution required in situations where KHI may 
become important. 

We also find that to capture strong shocks, like that 
expected in supernova explosions for example, a relatively 
high AV parameter, a AV , is required, which may be higher 
than what was been used previously for conventional galaxy 
formation simulations. We have applied this higher a AV to 
the new code, demonstrating that the new code is still able 
to model KHI, while simultaneously capturing the strong 
shock. 

It is also intriguing that our new code leads to similar 
results found when employing AMR codes to simulate the 
formation of a galaxy cluster within a cosmological context. 
We believe that the new AV and AC schemes significantly 
help to improve conventional SPH codes, and the resulting 
code should provide a means to simulate more accurately 
the formation and evolution of galaxies. 

In a forthcoming paper, we will carry out more realis- 
tic simulations of galaxy formation and evolution, including 
radiative cooling, star formation, SNe feedback and chemi- 
cal evolution, comparing and contrasting its behaviour with 
conventional SPH approaches. 
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